function fv=LS1(y1,x,r0,a,b,runmax,dmy,tolv,y2,frac)
%********************************************************
%** This program is to execute least square estimation **
%** for LSxtreg.m                                      **
%** Created by Min-Hua Huang, Jan 25, 2012             **
%** Model: Least Squares, by Fixing at Minimum, with   **
%**        Constrained Optimization                    **
%** y1: inividual-level dependent variable             **
%** x: covariate matrix                                **
%** r0: initial parameter value                        **
%** a: lower bound                                     **
%** b: upper bound                                     **
%** runmax: number of maximum iterations               **
%** dmy: number of dummy variables                     **
%** tolv: tolerance value                              **
%** y2: county-level dependent variable                **
%** frac: rescaled parameter, eta                      **
%********************************************************
format long; warning off all;
% basic setups
hr=boundary(y2,frac); b=b-hr(1); a=a-hr(2);
m=size(x,2)-1; record=ones(1,4*m+13)*-999; 
% initial values
lamda=ones(2*m+4,1)*10; kk=1; run=1; stop=1;
%****************** Record the initial results ***************************
record(1,1)=0;
record(1,2:m+2)=r0';
record(1,m+3:2*m+3)=zeros(1,m+1);
record(1,2*m+4:4*m+7)=lamda';
record(1,4*m+8)=(y1-x*r0)'*(y1-x*r0);
record(1,4*m+9:4*m+13)=ones(1,5)*-999;
% maximum and miminum vector of covariates 
xmax=zeros(m+1,1); xmin=zeros(m+1,1);
   for i=1:1:m+1
       xmax(i)=max(x(:,i));
       xmin(i)=min(x(:,i));
   end    
xlimit=[xmax xmin];   
while stop>0 && run<runmax
% vectors of indicator variables
v=vsign(r0);
% constraint vector
ci=cpar(r0,a,b,x,xlimit,v,dmy);
% gradient vector
df=-2*(x'*y1)+2*(x'*(x*r0));
% first derivative of the constraint vector
ai=aipar(r0,x,xlimit,v,dmy); 
%****************** Step 1:check KKT conditions ********************
% KKT
KKT1=df+ai'*lamda; check1=sum(abs(KKT1)>tolv);
KKT2=ci; check2=sum(KKT2>tolv);
KKT3=lamda; check3=sum(KKT3<-tolv);
KKT4=lamda'*ci; check4=sum(abs(KKT4)>tolv);
check=check1+check2; %+check3+check4;
   if check==0
      stop=0; 
   end
% If KKT conditions are not satisfied, then iteration continues.
%****************** Step 2:compute dk and lamdaQP ********************
% Hessian matrix
L=2*(x'*x);   
% solving the QP problem  
d=quadprog(L,df,ai,-ci);  
lamdaQP=(ai')\(-df-L*d);
%****************** Step 3:set new x and lamda ***********************
r1=r0+d;
lamda=lamdaQP;
%****************** Record the updated results ***********************
record(kk+1,1)=kk;
record(kk+1,2:m+2)=r1';
record(kk+1,m+3:2*m+3)=d';
record(kk+1,2*m+4:4*m+7)=lamdaQP';
record(kk+1,4*m+8)=(y1-x*r1)'*(y1-x*r1);
record(kk+1,4*m+9)=record(kk+1,4*m+8)-record(kk,4*m+8);
record(kk+1,4*m+10:4*m+13)=[check1 check2 check3 check4];
%****************** Step 4 and 5: return to Step 1 *******************
r0=r1;
kk=kk+1;
run=run+1;
end
minL=10^20; mint=-999;
for i=1:1:kk-1
    if record(i,4*m+8)<minL && record(i+1,4*m+11)==0
       minL=record(i,4*m+8);
       mint=i; 
    end
end
if mint~=-999 
   fv=[record(mint,2:m+2) record(mint+1,4*m+11) record(mint,4*m+8)];
else    
   fv=ones(1,m+3)*-999; 
end    
  